#d ;
cap log close ;
cap program drop nocoint ;
cap postclose nocoint ;
set more off ;
clear ;
set seed 1234567891 ;
cap mkdir "Set directory here" ;
global nocoint "Set directory here" ;
cap mkdir "$nocoint\figures" ;
cap mkdir "$nocoint\results" ;
cap mkdir "$nocoint\results\graphs" ;

program define nocoint, rclass ;
syntax [, i(real 0) j(real 0) k(real 0) l(real 0) m(real 0) n(real 0) p(real 1) q(real 1) s(real 0) obs(real 200)] ;
	#d ;
	drop _all ;
	set obs `obs' ;
	#d ;
	gen t = _n ; 
	tsset t ;
	matrix C = (1, 0 \  0, 1) ;
	drawnorm e u, cov(C) ;
	gen x=e if t==1 ;
	gen y=u if t==1 ;
	replace x= `q'*l.x + e if t>1 ; /*`q': memory of x.*/
	replace y= `p'*l.y + u if t>1 ; /*`p': memory of y.*/
	drop if t<=100 ; /*Discard the first 100 observations to eliminate the influence of initial values.*/
	
		/*Identify appropriate lags for the Dicky-Fuller test (y). */ 
		reg d.y l.y ; 
		predict resid_y, r ;
		wntestq resid_y ;
		local i=r(p) ; 
		drop resid_y ;
		while `i'<0.05 { ; /*The "while" command allows us to keep running the set of stata commands 
		                   enclosed in the braces until the p-value from the White noise Q test becomes 
		                   greater than 0.05. */
		local j =`j'+ 1 ; /*If the p-value is lower than 0.05, which is a sign of the presence of serial correlation, 
							`j' increases by one unit such that one more lag of ld.y is added in the equation below.*/
		reg d.y l.y l(1/`j')d.y ; 
		predict resid_y, r ;
		wntestq resid_y ;
		local i=r(p) ; 
		drop resid_y ;
		} ;
		
		/*Identify appropriate lags for the Dicky-Fuller test (x). */ 
		reg d.x l.x ; 
		predict resid_x, r ;
		wntestq resid_x ;
		local k=r(p) ; 
		drop resid_x ;
		while `k'<0.05 { ; 
		local l=`l'+ 1 ; 
		reg d.x l.x l(1/`l')d.x ; 
		predict resid_x, r ;
		wntestq resid_x ;
		local k=r(p) ; 
		drop resid_x ;
		} ;
		
		/*Identify appropriate lags for the Engle-Granger (EG) test. */
		reg y x ;
		predict resid_eg_1, r ;
		reg d.resid_eg_1 l.resid_eg_1, nocons ;
		predict resid_eg_2, r ;
		wntestq resid_eg_2 ;
		local m=r(p) ;
		drop resid_eg_1 resid_eg_2 ;
		while `m'<0.05 { ;
		local n = `n'+ 1 ; 
		reg y x ;
		predict resid_eg_1, r ;
		reg d.resid_eg_1 l.resid_eg_1 l(1/`n')d.resid_eg_1, nocons ;
		predict resid_eg_2, r ;
		wntestq resid_eg_2 ;
		local m=r(p) ;
		drop resid_eg_1 resid_eg_2 ;
		} ;
		

		/*Augmented Dicky-Fuller test */
		dfuller y, lags(`j') ; /*dfuller: dfuller=Dickey-Fuller (DF) Test
					H0: I(1) or y is unit root.
					H1: (p-1)<0 or y is stationary.
					The appropriate lags `j' is imputed so as to address 
					the issues of serial correlation.
					If p-value is greater than 0.05, it indicates that y
					has a unit root.
					Note that rejecting the null hypothesis does not mean
					we can accept H0 (Woodridge 2009, 633)!
					dfuller by default regresses d.y on l.y. 
					How many lags to be included often depends on 
					the frequency of the data.*/
		
		return scalar p_value_df_y=r(p) ; /*Save p-value to see if the 
										estimated coefficient on l.y is 
										statistically significant at the .05 level 
										given the Dicket-Fuller Distribution.*/
		
		dfuller x, lags(`l') ;
		return scalar p_value_df_x=r(p) ; 
		
		/*Engle-Granger (EG) Test */
		egranger y x, lags(`n') ; /*egranger: Engle-Granger (EG) Test
						H0: residual being nonstationary (meaning that y and x are NOT cointegrated)
						H1: residual being stationary (indicating cointegration)
						egranger by default implements the following analysis:
						. regress y x
						. predict double resid, res
						. regress D.resid L.resid, nocons 
						Again, as with the DF test we can increase lags of 
						resid to account for serial correlation.*/
		return scalar Zt= e(Zt) ; /*Save test t-stat for l.resid. */
		return scalar cv5 = e(cv5) ; /*Save the 5% critical value*/
		
			
		return scalar tstat_gr2 = abs(_b[L1._egresid]/_se[L1._egresid]) ; 
		
		
		/*Error Correction Model (ECM)*/
		reg d.y l.y d.x l.x ;
		return scalar beta_ly=_b[L.y] ; /*Save these to check bias in our esitmated coefficients.*/
		return scalar beta_dx=_b[D.x] ;
		return scalar beta_lx=_b[L.x] ;
		return scalar tstat_dx = abs(_b[D.x]/_se[D.x]) ; 
		return scalar tstat_lx = abs(_b[L.x]/_se[L.x]) ; 
		return scalar tstat_ly = _b[L.y]/_se[L.y] ; 
		matrix V=e(V) ;
		matlist V ;
		local LRM = _b[L.x]/abs(_b[L.y]) ;
		local var_LRM = (1/_b[L.y]^2)*V[3,3]+(_b[L.x]^2/_b[L.y]^4)*V[1,1]
		-2*(_b[L.x]/_b[L.y]^3)*V[3,1] ;
		local se_LRM=sqrt(`var_LRM') ;
		return scalar tstat_LRM = abs(`LRM'/`se_LRM') ;
		return scalar p=`p' ;
		return scalar q=`q' ;
		return scalar obs=`obs' ;
		
end ;
nocoint ;

#d ;
postfile nocoint  obs p q 
				fsig_df_y fsig_df_x 
				fsig_gr_mac_c fsig_gr_t fsig_gr_t_2 
				fsig_ecm_mac_nc fsig_ecm_mac_c fsig_ecm_t fsig_ecm_t_2 
				fsig_lx_0 fsig_lx_0_df_y fsig_dx_0 fsig_dx_0_df_y 
				fsig_LRM_0 fsig_LRM_0_df_y 
				fsig_lx_ecm_mac_nc fsig_LRM_ecm_mac_nc 
				fsig_lx_ecm_mac_c fsig_LRM_ecm_mac_c 
				fsig_lx_ecm_mac_nc_df_y fsig_LRM_ecm_mac_nc_df_y 
				fsig_lx_ecm_mac_c_df_y fsig_LRM_ecm_mac_c_df_y 
				fsig_lx_ecm_t fsig_LRM_ecm_t 
				fsig_lx_ecm_t_2 fsig_LRM_ecm_t_2 
				fsig_lx_ecm_t_2_dy fsig_LRM_ecm_t_2_dy 
				fsig_lx_gr_mac_c fsig_LRM_gr_mac_c 
				fsig_lx_gr_mac_c_df_y fsig_LRM_gr_mac_c_df_y 
				fsig_lx_gr_t fsig_LRM_gr_t 
				fsig_lx_gr_t_2 fsig_LRM_gr_t_2 
				beta_mean_ly beta_mean_dx beta_mean_lx 
using "$nocoint\results\results_nocoint", replace ; 
foreach obs of numlist 
130 160 200 { ;
foreach q of numlist 
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.99 1 { ;
foreach p of numlist 
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.99 1 { ;
				simul p=r(p) q=r(q) obs=r(obs) 
				p_value_df_y=r(p_value_df_y) p_value_df_x=r(p_value_df_x) 
				Zt=r(Zt) cv5=r(cv5) 
				tstat_dx=r(tstat_dx) tstat_lx=r(tstat_lx) tstat_ly=r(tstat_ly) tstat_LRM=r(tstat_LRM) 
				beta_ly=r(beta_ly) beta_dx=r(beta_dx) beta_lx=r(beta_lx), 
				reps(2000): 
				nocoint, obs(`obs') p(`p') q(`q') ;
				local obs_n =`obs'-100 ;
				local obs_df =`obs'-100-4 ; /*Four explanatory variables including a constant.*/
				local obs_ecm_mac_nc =`obs'-100-3 ; /*Three explanatory variables excluding a constant.*/
				local cv_t_one = invttail(`obs_df',0.05) ;
				local cv_t_two = invttail(`obs_df',0.025) ;
				local cv_gr_mac_c = -3.33613-6.1101/`obs_n'-6.823/(`obs_n'^2) ; /*From MacKinnon (2010)*/
				local cv_ecm_mac_c = -3.2145-3.21/`obs_df'-2/(`obs_df'^2) +17/(`obs_df'^3) ; /*From MacKinnon and Ericsson (2002)*/
				local cv_ecm_mac_nc = -2.5937-1.53/`obs_ecm_mac_nc'-0.8/(`obs_ecm_mac_nc'^2)-24/(`obs_ecm_mac_nc'^3) ; /*From MacKinnon and Ericsson (2002)*/
				
				/*Results from the DF test. */
				sum p_value_df_y if p_value_df_y>0.05 ; 
				local fsig_df_y = r(N)/20 ; /*The total number of rejecting H0 of I(1) from the DF test.*/
				sum p_value_df_x if p_value_df_x>0.05 ; 
				local fsig_df_x = r(N)/20 ;
				
				/*Examne how often we find evidence of spurious regression for all series. */
				sum tstat_lx if tstat_lx>=`cv_t_two' ; 
				local fsig_lx_0 = r(N)/20 ;
				sum tstat_dx if tstat_dx>=`cv_t_two' ; 
				local fsig_dx_0 = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' ; 
				local fsig_LRM_0 = r(N)/20 ;
				
				sum tstat_lx if tstat_lx>=`cv_t_two' & p_value_df_y>0.05 ; 
				local fsig_lx_0_df_y = r(N)/20 ;
				sum tstat_dx if tstat_dx>=`cv_t_two' & p_value_df_y>0.05 ; 
				local fsig_dx_0_df_y = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' & p_value_df_y>0.05 ; 
				local fsig_LRM_0_df_y = r(N)/20 ;
				
				/*ECM test for cointegration*/
				sum tstat_ly if tstat_ly<=`cv_ecm_mac_nc' ; 
				local fsig_ecm_mac_nc = r(N)/20 ;
				sum tstat_ly if tstat_ly<=`cv_ecm_mac_c' ;
				local fsig_ecm_mac_c = r(N)/20 ;
				sum tstat_ly if tstat_ly<=-`cv_t_one' ; 
				local fsig_ecm_t = r(N)/20 ;
				sum tstat_ly if tstat_ly<=-`cv_t_two' ; 
				local fsig_ecm_t_2 = r(N)/20 ;
				
				/*Examine how often we find evidence of spurious regression, only for the series for which we find evidence of cointegration by the ECM test. */
				sum tstat_lx if tstat_lx>=`cv_t_two' & tstat_ly<=`cv_ecm_mac_nc' /*MacKinnon critical value without a constant*/ ; 
				local fsig_lx_ecm_mac_nc = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' & tstat_ly<=`cv_ecm_mac_nc'; 
				local fsig_LRM_ecm_mac_nc = r(N)/20 ;
				
				sum tstat_lx if tstat_lx>=`cv_t_two' & tstat_ly<=`cv_ecm_mac_c' /*MacKinnon critical value with a constant*/ ; 
				local fsig_lx_ecm_mac_c = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' & tstat_ly<=`cv_ecm_mac_c' ; 
				local fsig_LRM_ecm_mac_c = r(N)/20 ;
				
				sum tstat_lx if tstat_lx>=`cv_t_two' & tstat_ly<=-`cv_t_one' ; 
				local fsig_lx_ecm_t = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' & tstat_ly<=-`cv_t_one' ; 
				local fsig_LRM_ecm_t = r(N)/20 ;
				
				sum tstat_lx if tstat_lx>=`cv_t_two' & tstat_ly<= -`cv_t_two' ; 
				local fsig_lx_ecm_t_2 = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' & tstat_ly<= -`cv_t_two' ;
				local fsig_LRM_ecm_t_2 = r(N)/20 ;
				sum tstat_lx if tstat_lx>=`cv_t_two' & tstat_ly<= -`cv_t_two' & p_value_df_y>0.05 ; 
				local fsig_lx_ecm_t_2_dy = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' & tstat_ly<= -`cv_t_two' & p_value_df_y>0.05 ;
				local fsig_LRM_ecm_t_2_dy = r(N)/20 ;
				
				sum tstat_lx if tstat_lx>=`cv_t_two' & tstat_ly<= `cv_ecm_mac_nc' & p_value_df_y>0.05  ; 
				local fsig_lx_ecm_mac_nc_df_y = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' & tstat_ly<= `cv_ecm_mac_nc' & p_value_df_y>0.05  ;
				local fsig_LRM_ecm_mac_nc_df_y = r(N)/20 ;
				
				sum tstat_lx if tstat_lx>=`cv_t_two' & tstat_ly<= `cv_ecm_mac_c' & p_value_df_y>0.05  ; 
				local fsig_lx_ecm_mac_c_df_y = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' & tstat_ly<= `cv_ecm_mac_c' & p_value_df_y>0.05   ;
				local fsig_LRM_ecm_mac_c_df_y = r(N)/20 ;
				
				
				/*EG test for cointegration */
				sum Zt if Zt<=`cv_gr_mac_c' ; 
				local fsig_gr_mac_c = r(N)/20 ; /*The total number of rejecting H0 of nonstationarity from the EG test, thus detecting evidence of cointegration.*/
				sum Zt if Zt<=-`cv_t_one' ;
				local fsig_gr_t = r(N)/20 ; 
				sum Zt if Zt<=-`cv_t_two' ;
				local fsig_gr_t_2 = r(N)/20 ;
				
				/*Examine how often we find evidence of spurious regression, only for the series for which we find evidence of cointegration by the EG test. */
				sum tstat_lx if tstat_lx>=`cv_t_two' & Zt<=`cv_gr_mac_c' /*MacKinnon critical value*/ ; 
				local fsig_lx_gr_mac_c = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' & Zt<=`cv_gr_mac_c' ; 
				local fsig_LRM_gr_mac_c = r(N)/20 ;
				
				sum tstat_lx if tstat_lx>=`cv_t_two' & Zt<=-`cv_t_one' ; 
				local fsig_lx_gr_t = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' & Zt<=-`cv_t_one' ; 
				local fsig_LRM_gr_t = r(N)/20 ;
				
				sum tstat_lx if tstat_lx>=`cv_t_two' & Zt<= -`cv_t_two' ; 
				local fsig_lx_gr_t_2 = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' & Zt<= -`cv_t_two' ; 
				local fsig_LRM_gr_t_2 = r(N)/20 ;
				
				sum tstat_lx if tstat_lx>=`cv_t_two' & Zt<=`cv_gr_mac_c'  & p_value_df_y>0.05   ; 
				local fsig_lx_gr_mac_c_df_y = r(N)/20 ;
				sum tstat_LRM if tstat_LRM>=`cv_t_two' & Zt<=`cv_gr_mac_c' & p_value_df_y>0.05   ; 
				local fsig_LRM_gr_mac_c_df_y = r(N)/20 ;
				
				/*Save coefficients*/
				sum beta_ly ; 
				local beta_mean_ly = round(r(mean), .001) ; 
				
				sum beta_dx ; 
				local beta_mean_dx = round(r(mean), .001) ; 
				
				sum beta_lx ; 
				local beta_mean_lx = round(r(mean), .001) ; 
				
				
				post nocoint (`obs') (`p') (`q') 
				(`fsig_df_y') (`fsig_df_x') 
				(`fsig_gr_mac_c') (`fsig_gr_t') (`fsig_gr_t_2') 
				(`fsig_ecm_mac_nc') (`fsig_ecm_mac_c') (`fsig_ecm_t') (`fsig_ecm_t_2') 
				(`fsig_lx_0') (`fsig_lx_0_df_y') (`fsig_dx_0') (`fsig_dx_0_df_y') 
				(`fsig_LRM_0') (`fsig_LRM_0_df_y') 
				(`fsig_lx_ecm_mac_nc') (`fsig_LRM_ecm_mac_nc') 
				(`fsig_lx_ecm_mac_c') (`fsig_LRM_ecm_mac_c') 
				(`fsig_lx_ecm_mac_nc_df_y') (`fsig_LRM_ecm_mac_nc_df_y') 
				(`fsig_lx_ecm_mac_c_df_y') (`fsig_LRM_ecm_mac_c_df_y') 
				(`fsig_lx_ecm_t') (`fsig_LRM_ecm_t') 
				(`fsig_lx_ecm_t_2') (`fsig_LRM_ecm_t_2') 
				(`fsig_lx_ecm_t_2_dy') (`fsig_LRM_ecm_t_2_dy') 
				(`fsig_lx_gr_mac_c') (`fsig_LRM_gr_mac_c') 
				(`fsig_lx_gr_mac_c_df_y') (`fsig_LRM_gr_mac_c_df_y') 
				(`fsig_lx_gr_t') (`fsig_LRM_gr_t') 
				(`fsig_lx_gr_t_2') (`fsig_LRM_gr_t_2') 
				(`beta_mean_ly') (`beta_mean_dx') (`beta_mean_lx') ;
				} ;
				} ;
				} ;
				postclose nocoint ;	

/*Replicate Table A-1*/
#d ;
use "$nocoint\results\results_nocoint", clear ;
keep obs p q fsig_lx_0 fsig_lx_ecm_t_2_dy fsig_lx_ecm_mac_c_df_y ;
order obs p q fsig_lx_0 fsig_lx_ecm_t_2_dy fsig_lx_ecm_mac_c_df_y ;
keep if q>0.85 ;
keep if p>0.85 ;
local n2=`n'-100 ;
outsheet using "$nocoint\results\Table_A-1.csv", comma replace ;

/*Replicate Table A-2*/
#d ;
use "$nocoint\results\results_nocoint", clear ;
keep obs p q fsig_LRM_0 fsig_LRM_ecm_t_2_dy fsig_LRM_ecm_mac_c_df_y ;
order obs p q fsig_LRM_0 fsig_LRM_ecm_t_2_dy fsig_LRM_ecm_mac_c_df_y ;
preserve ;
keep if q>0.85 ;
keep if p>0.85 ;
local n2=`n'-100 ;
outsheet using "$nocoint\results\Table_A-2.csv", comma replace ;

